Fully synchronous solutions and the synchronization phase 
transition for the finite-TV Kuramoto model 
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Abstract 

We present a detailed analysis of the stability of synchronized solutions to the Kuramoto system of 
oscillators. We derive an analytical expression counting the dimension of the unstable manifold associated 
to a given stationary solution. From this we are able to derive a number of consequences, including: 
analytic expressions for the first and last frequency vectors to synchronize, upper and lower bounds on 
the probability that a randomly chosen frequency vector will synchronize, and very sharp results on the 
large N limit of this model. One of the surprises in this calculation is that for frequencies that are Gaussian 
distributed the correct scaling for full synchrony is not the one commonly studied in the literature — rather, 
there is a logarithmic correction to the scaling which is related to the extremal value statistics of the random 
frequency vector. 

1 Introduction 

1.1 History of Kuramoto model 

The study of synchronization of coupled nonlinear oscillators has a history that spans several centuries, 
starting with Huygens' observation of synchronizing pendulum clocks |3 12]. There has been a great body 
of work throughout this history studying such synchronization phenomena; for reviews see pl|[25)[30) . 
There are a wide variety of such models derived from a diverse collection of mathematical and scientific 
questions, including pulse-coupled ] 13 ,20] and conservative |7,8| models. These models arise in a number 
of contexts, including numerous biological applications 114JI17L and exhibit a diversity of behaviors. 

A fundamental model of synchronization, however, is the case where we consider independent oscillators 
connected through dissipative coupling, and the fundamental question is |25 26|: how do independent 



oscillators with different frequencies adjust themselves to produce a collective mode? 
The system of coupled ordinary differential equations 



dt 



3 



Sill 
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was first proposed by Kuramoto ] 15 161 as a model for the synchronization of oscillators; this model, and 



variants, are widely known as the Kuramoto model. Since this time, the Kuramoto model has been a fun- 
damental model for many types of systems exhibiting synchronization |2]|6]|11J|27| and related phenomena 
such as flocking |[9p0|. For reviews, see II 26]. Most work in this area has analyzed the Kuramoto model in 
the continuum limit, where the number is oscillators N is formally allowed to go to infinity and the sum- 
mation is replaced by an appropriate integral. Such formal calculations provide a tremendous amount of 
physical insight into the problem; however, it has proven to be difficult to make these approaches rigorous. 

An alternative approach, attributed by Strogatz |26) to a series of lectures by Kopell, is to analyze the 
finite- N problem carefully and then take the A" — > oo limit. There have been a few papers which have 
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established rigorous results on the existence and stability of synchronized solutions ] 18,. 19 28 , 29 J but to 
date this the finite- N problem has been only partially understood. In this paper we carry out a substantial 
portion of this program. We focus on the case of full synchrony where all of the oscillators rotate with the 
same angular frequency. We derive new characterizations of the stable regions, which permit a detailed 
understanding of the region in frequency space where full synchronization occurs. In particular we are 
able to identify a large subset of the stable region, in the form of the Voronoi cell of a well-understood 
high dimensional lattice. Using this we find upper and lower bounds on the probability that the system 
undergoes full synchronization. 



1.2 Problem Formulation 

The standard formation of the Kuramoto problem is as follows: Consider a weighted (directed) graph 
G = (V, E) and denote 7^ > as the weight of edge i j. For any lj € R N , define the dynamical system 
on 9 eT N by 

^9, ^uj. + Y^ Hi sin^ - (1.1) 

3 

One of the more common choices of interaction graph is the symmetric all-to-all graph where we assume 
that all the 7^ are equal. We thus consider the system 

^ i =cj i + 7 ^sin(^-^), (1.2) 

3 

or, defining the function f : T" — > R n coordinatewise as 

fi{O)=Y,M0i-0i)> (1-3) 



^0 = « + 7f(0). (1.4) 



we can write 1 1.2 1 as 



Remark 1.1 Note that we have not rescaled the coupling coefficient in \1.2) ; the typical scaling chosen for Kuramoto 
is 7 /AT. We will refer to this scaling below as the "classical scaling". One of the results of this paper is that the classical 
scaling is only the correct scaling for certain problems, and not for others. For the problem where the frequencies Wj 
are chosen from a Gaussian distribution, for instance, we shall see that the correct scaling differs from the classical one 
by a logarithmic term. 

Remark 1.2 It is worth noting that this flow is a gradient flow, if one considers the angles as lying in the covering 
space rather than T N . The flow can be written as 

— = VL 

dt 

where the energy L is given by 

L= (w, 0) + 7 ^ cos (° l ~ 6 i)~ N - 

hi 

Here the constant —N is chosen for convenience and obviously doesn't influence the dynamics. It is common to work 
with the "order parameter" R 



cosie, -6j), 

3 



2 



which gives the energy function as 

L = (lj,8)+R 2 (6). 

Thus the Kuramotoflow tries to maximize an energy given by a sum of two terms. The first term acts to align the flow 
with the frequency vector oj, while the second acts to increase the order parameter. All of the structure that arises in 
Kuramoto model is due to a competition between these two effects. 



The fundamental question we consider is the question of whether 1 1.2 1 admits a fully synchronous 
solution, as we now define: 

Definition 1.1 For a given value of the frequency vector u we say that the Kuramoto model exhibits full synchrony 
if \1.2) with 7 = 1 admits a stationary solution which is asymptotically stable, i.e. we have a solution to the 
equation 

f (0) = -u (1.5) 

such that, if we define the jacobian matrix 

J = Vof, (1-6) 

then J is negative semi-definite with a one-dimensional kernel. If such a solution exists, we call it a fully synchronous 
solution. 

Remark 1.1 It is not hard to see that, for fully synchronous solutions, we will have 

max sup \\0i(t) — 6 j(t)\\ < oo. 

i'i te[o,oo) 



However, note that due to the continuous symmetry of \1.2) , these fully synchronous solutions can be moving in 
the same co-rotating frame. The Kuramoto model also permits partially synchronous solutions where a large mass of 
the oscillators are fixed relative to each other, and some subset precesses relative to them. We do not consider such 
solutions in this work, q.v. the discussion in Section^below. 

The results of the paper are organized as follows. In Section|2]we give an almost complete characteriza- 
tion of the u> which give rise to fully synchronous solutions; in particular, we show that this set is convex, 
identify many important points on the boundary of this set, and explicitly determine the convex hull of 
these points. We would like to characterize "how easy" it is for a Kuramoto problem to fully synchronize, 
but there are many senses in which one can pose this problem; we consider two different cases below. 

For example, we could pose the following question: given u> € R , define 7*(w) as the minimal 7 for 
which | |1.4[ | has a fully synchronous solution. Of course the presence of synchrony is invariant under a 
rescaling in time, so one may as well assume that = 1. The motivates the definition 

Definition 1.2 Define the lower and upper critical couplings 7 m i„ (N) and 7 max (./V) in the following way: 

7min (A0 := inf 7 *(w), 7 max (iV) := sup 7 *(w) (1.7) 



Roughly 7 m i n (iV) characterizes the minimal coupling constant we should choose to get full synchrony, and 
the values of u> for which the infimum is attained (by symmetry there are many) can be thought of as the 
"most easily synchronizable frequency(s)"; analogously 7 max (A r ) characterizes the most difficult frequen- 
cies, to synchronize, and gives the largest necessary coupling constant. For coupling constants above this 
value all frequency vectors synchronize. The quantities 7 m i n and 7 max are closely related to the quantity K c 
studied by Verwoerd and Mason in their work |28||29) : 7 m i n (iV) is essentially the minimum of K c over the 
unit circle, and 7 max the maximum. 

We study 7min(A0, 7max(A0 in Section|3| and give bounds which show the different asymptotic scalings 
of these quantities. We will show that in the limit N — > 00, 7 m i n (^V) ~ N~ 3 / 2 and 7 max (^V) ~ AT -1 . In 
particular, synchrony "turns on" much earlier than the classical scaling, although it "fills up" in a way 
consistent with the classical scaling. 
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A different method of choosing u), which is very common in the literature, is to assume that the entries 
of u> are chosen independently from some probability distribution, i.e. that the uj n are iid random variables. 
In Section |4] we consider this question; stated precisely, we proceed as follows. Fix 7 > and N, and 
choose Lu n independently from a Gaussian distribution with mean zero and unit variance, i.e. the u> n are 
independent N(0, 1). Define V S ync(l,N) as the probability of (172} havi ng a fully synchronous solution. 
Define ip(N) := ^2 \n{N)/(N + 1), and what we show below in Theorem 



4.1 



is 




lim V sync (6ip(N),N) = ^' r J (1.8) 

TV— »oo 

In particular, this implies that V sy nc{l/N, N) — > for any 7, so that in the classical scaling the probability of 



full synchrony is zero. In fact, we will show in Proposition 4.5 that the probability of full synchrony decays 
to zero exponentially fast as N —> 00. We show that this anomalous scaling is closely related to the extreme 
value statistics for a Gaussian distribution, and is to be expected for any distribution of frequencies which 
is not of compact support. 

Finally, in Section [6] we finish with some comments about the different choices of scaling and why they 
arise and present a few open questions. 

2 Characterization of the Stable Set 



In this section we write down a relatively complete description of the set of u> for which 1 1.2 1 admits a stable 
fully synchronous solution. We do this by proving an index theorem which counts the dimension of the 
unstable manifold to any stationary solution. The stable synchronous solutions are then those for which 
the unstable manifold is zero dimensional. 

2.1 Notation 

We consider the finite N Kuramoto model with uniform sine coupling on the complete graph, which we 
repeat 

d 



-O = w + 7f(0) (L4I 



Note that 

N 



for all 9, due to a telescoping sum. Therefore we have 

d 



dt 



^0 4 = 5> 4 =:a (2.1) 



This means that, if f2 7^ 0, the center of mass of the system precesses around the circle with a constant 
velocity. Using the change of variables 6>; = B L — N~ x £Lt puts us into a corotating frame and allows us to 
assume without loss of generality that Y^i^i = 0- 

The function f : T" — > R n is a natural map from the configuration space T n to the frequency space R™. 
Assuming some convexity conditions which will be established later, the frequency vector uj is the convex 
dual variable to the angle vector 6. It will also be convenient to think of this as a map f : T n_1 — > W 1 ^ 1 = 
M n /{1, 1, . . . , 1} — this can be done, for example, by fixing a single 9 n . 

The Jacobian of this mapping controls the stability of the synchronized state, and it is crucial to have 
a good understanding of set of parameter values for which the Jacobian is negative semi-definite. This 
motivates the following definition: 
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Definition 2.1 We define Sg to be the set of configurations for which the Jacobian || is negative semi-definite with 
a one dimensional kernel. The boundary of this set is obviously the set of configurations for which the Jacobian || is 
negative semi-definite with a kernel of dimension two or more. 

Throughout the paper calligraphic capital letters will denote sets of particular interest, with a subscript 
of 8 denoting sets in the configuration space T" -1 and a subscript u denoting the corresponding sets in the 
frequency space. For instance if Sg denotes the set of asymptotically stable stationary configurations then 
S^, the set of frequencies admitting an asymptotically stable stationary configuration, is simply the image 
of Sg under the map f: 

f : Sg n- 

It is clear that S^ is the important object for studying synchronization: all questions about the probability 
of (full) synchrony are questions about the size of S^, in some measure. One of the key ingredients in this is 
a good characterization of S^, . 



2.2 Index Theorem 

In this section we prove an index theorem which counts the dimension of the unstable manifold to the 
synchronized state. The basic idea of this section is that the Jacobian takes a relatively simple form, as it 
can be written as a rank two perturbation of a diagonal matrix. An application of a rank-two perturbation 
formula gives a straightforward count of the number of positive eigenvalues. 

Definition 2.2 Given a matrix A, we define the three quantities n + (A), ri-(A), no(A) as the number of eigenvalues 
of A with positive, negative, and zero real parts, respectively. We will refer to these as the indices of A. (Given a 
vector field f with fixed point 0* , we will abuse notation and refer to the indices of 6* when we mean the indices of the 
jacobian of the vector field at 0.) 

A straightforward calculation gives the following expression for the Jacobian matrix 



Me) 



dfi 
80: 



(0) 



cos 



Using the cosine angle addition formula, we have 

J = -£> + v<g)v- 



where the vectors v and w are defined by 

sin(0 2 ) 



cos(# fe 



i), i = j- 



\ sin(0„) / 



D is the diagonal matrix 



(2.2) 



(2.3) 



Di 



E' 



/ cob(0i) 
cos(6> 2 ) 

V cos(6»„) J 
h), 



(2.4) 



(2.5) 



and £§) denotes the usual outer product (a ® b)^ = a,6j. 

The Jacobian is a positive semi-definite rank two perturbation of a diagonal matrix, a fact which figured 
in the analysis of p9) . The spectrum of a low rank perturbation of an operator can be computed explicitly in 
terms of the spectrum of the unperturbed operator together with some information about the inverse of the 
(unperturbed) operator, a fact usually known as the Aronszajn-Krein formula |24|. Here the unpertubed 
operator is diagonal and can thus be trivially inverted, and one can get a reasonably complete description 
of the spectrum. 
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Thus we proceed as follows. Since the full Jacobian is a rank two perturbation of the diagonal, the 
indices of J and D can differ by at most two. We will first compute the index of D, then compute the 
difference of indices explicitly. 



Definition 2.3 Given a configuration 9, we define the complex order parameter 

N 

N 



1 - 

fe^-^e 19 ". (2.6) 



n=l 



Proposition 2.1 Using the 0(1) invariance, rotate the configuration so that the order parameter is on the positive 
real axis. Then we have 

n+{-D) = #{6 i \cos{e i ) < 0} 

Proof. Since D is diagonal we just have to count the number of negative entries on the diagonal. We 
have 

Da = cos (^ _ B *) = cos( - 9 ^ cos{6 k ) + sin(0 fc ) ^ sin(^) (2.7) 

k k 

which we recognize as the (1,0) component of the order parameter rotated through angle 0j. The number 
of times this is negative is the number of angles in the range (tt/2, 3tt/2). ■ 

Remark 2.1 Proposition 2.1 is equivalent to the calculations done in Section 4 of [19]. 

Since the perturbation is rank two, and one of the eigenvalues of the Jacobian is zero, if D is negative- 
definite then there is only one eigenvalue whose sign is to be determined. 

Theorem 2.2 Suppose D is invertible. Define the following quantity: 
Then we have 

n + (J)=n + (-D) + {°' T<2 ' (2.9) 

[1, T > 2. 

Proof. We proceed by a homotopy argument in the spirit of the Birman-Schwinger Principle (see, for 
instance, the discussion on page 98 of Volume 4 of the text of Reed and Simon p2|). Consider the one 
parameter family of operators 

J n = -D + Tj(v ® v + w ® w) (2.10) 

Clearly J = —D and J\ = J. Since v<g)v + w<g>wis non-negative, all eigenvalues are non-decreasing func- 
tions of r). We can detect eigenvalues crossing from the left half-plane to the right half-plane by detecting 
changes in the dimension of the kernel of J v . has non- trivial kernel if there is a vector x with J, ; x = 0, or 

— Z?x + ry((v, x) v + (w, x) w) = 0. 

Solving for x gives 

x = ?7((v,x)D" 1 v+(w,x) J D" 1 w). (2.11) 
Taking the inner product with v and w gives 

(v, x) = 77«v, x) (v, D-\) + (w, x> (v, D^w)), (2.12) 
(w, x) = n((v, x) (w, D- X v) + (w, x) (w, D _1 w)). (2.13) 
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This can be thought of as a linear system for the quantities 



ai = (v, x) , a 2 = (w, x) , (2.14) 

which can be written in matrix form as 

Mtl {a 2 )-{ r,(w,D-iv) r, (w, D^w) - 1 ){ a 2 ) ~ { J ' (Z15) 



There exists a nontrivial solution to ( 2.15[ (i.e., M v has a nontrivial kernel) if and only if J r) itself has a 



nontrivial kernel. Thus we want to count zeroes of det(M n ) for r\ € (0, 1). Initially Mq = —I 2X 2 has two 
negative eigenvalues. At r\ = 1 we have that Mi is singular so that one of the eigenvalues is zero. M is a 
linear matrix value function of r\, and the linear term 

(w,D _1 v) (w^-D^w) 

is positive definite, so the eigenvalues of M are increasing functions of 77. This implies that all eigenvalue 
crossings are transverse left to right and so the number of positive eigenvalues of the full matrix is deter- 
mined by the non-zero eigenvalue of M\. If this is negative there have been no eigenvalue crossings, and if 
this is positive there has been one. Since one eigenvalue is zero the other is equal to the trace, which is 

Tr(Mx) = -1 + (v, D _1 v) - 1 + (w, D _1 w) = (l, D^l) - 2 = r - 2. 

Thus if this quantity is positive a single eigenvalue has crossed into the right half -plane and n+(— D) = 
n + (J) + l. If this quantity is negative then no eigenvalues have crossed and the count is the same, n + (—D) = 
n+(J). 



Remark 2.1 Since the Jacobian matrix for the Kuramoto model is of the form of a graph Laplacian, the Kir choff matrix 
tree theorem can be applied. This expresses the product of the non-zero eigenvalues in terms of a sum over spanning 
trees of a product over edges in the spanning tree of the edge weights. In our case the edge weights may not be positive, 
but the standard proofs of the matrix tree theorem still hold in this case. Our formula above suggests that the product 
over the non-zero eigenvalues of the Jacobian should be proportional to the quantity 

Tr(Afi) = r - 2. 

One can check that this is the case, and actually compute the constant of proportionality. This gives the following 
unusual combinatorial identity 

^ ll cos ^) = y C0S (e- ~e ) ■ (Z16) 

Here S is the set of all spanning trees on N points (of which there are N N ~ 2 ), T is an element ofS, a spanning tree, 
e is an edge in the spanning tree, and if the edge e connects vertices I and m then 6 e = 61 — 6 m . The left-hand side 
of \2.16\ comes from the matrix tree theorem, and the righthand side comes from directly evaluating the product of the 
non-zero eigenvalues using arguments similar to those given above. This identity is not used in the current paper - 



we really only need the formula derived in Theorem 2.2 — but the existence of such a formula is extremely interesting. 



It suggests that there may be much more algebraic structure in this problem than is readily apparent. 

Having a concise index count for the Jacobian makes it possible to begin to make an analytical descrip- 
tion of this region. Our first observation is: 

Lemma 2.3 Let Se denote the region in configuration space where the there exists a stable synchronized solution. 
This region is simply connected and is given by the connected component o/t _1 ([1, 2)) that contains the origin. 
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Proof. To see that the stability region is simply connected, and that it is the component containing the 
origin, we give an explicit stability-preserving retraction of any stable steady state onto the stable steady 
state 6 = 0. Consider, for example, the retraction 9(() = £0. When £ = 1, we have the original steady state 
6, and when £ = 0, we have the stable steady state 9, = 0. 

The derivative of the Jacobian with respect to the homotopy parameter is given by 

dJ\ f (9 J -9 i )sm(((9 3 -9 l )) i+3 
dOa I -E fc (0fe-0i)sm(C(0fc-0i)). i=j 

A necessary condition for stability is that \9. L — 0j\ < tt for all i, j. This implies that (9j — 9i) sin(£(6>j — 9i)) 
is a positive quantity thus that is positive definite. Thus as the homotopy parameter Q is decreased 
from 1 to the eigenvalues of the Jacobian decrease, and all stable stationary states are contractible onto the 
state 9 = in a way that increases stability. This establishes the simple connectedness of Se and thus, by 
continuity, S^. 



We need something stronger than Lemma 2.3 In order to derive lower bounds on quantities related 
to the probability of stabilization we would like a much stronger result: namely that the stable region is a 
convex set. This is the content of the next few results. 

Definition 2.4 Define Ki(0) = J2j cos(6j — 9i). Define Se as the subset ofT n for which we have k, ; > and 
t = Y17=i K 7 1 — ^' an d S<*> as the image ofSe under the mwp f . 



Remark 2.2 Note that by Proposition 2.1 and Theorem 2.2 is precisely the set of oj such that uj = {(9) for some 



0, and the Jacobian \2.3\ is negative semi-definite; in short, those uifor which we have fully synchronous solutions. 

The above characterization of the stable region is difficult to use in practice, since one needs to find the 
image under a reasonably complicated map of a region defined by a transcendental equation. The following 
lemma shows that the boundary of this region is actually a connected piece of an algebraic variety. 

Lemma 2.4 The stable frequency set S w can be defined algebraically as follows: It is the solution set of 

K 2 i+ u? = ^ Kj (2.17) 
i 

n t > (2.18) 



E 1 



— < 2 (2.19) 

Ki 

5> i= (2.20) 
Since all ^ > 0, ( 2.19\ implies n t > \for all i. Also note that summing the first equation over i shows 

k 2 +u> 2 = N(1,k) 

or, equivalently 

1 1 4 

3 

and thus the entire stable region lies in a sphere of radius uj < 

Proof. The main observation is the following: squaring Ki and cjj and adding them gives 

• 2 1 ■ - 2 - cos(0j - 9,) cos(0,-/ - 6i) + sin(6j - 9. L ) sin(^- ° A 



"3 

3,3' 



Y, COS ( 6 3 

3,3' 
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or 



Fj(K,w) 



3=1 



Computing the Jacobian of the function F we find that 



dF 



I 2ki - 1 

-1 



-1 

2n 2 - 1 



-1 



-1 
-1 

-1 



0, for all i. 



(2.21) 



(2.22) 



2«n-l / 



It is straightforward to check that the determinant of this Jacobian is 



det( 



9k' 



yn—1 



N 

n< 

i=l 



E 1 



All of the principal minors are of the same form, so it is easy to check that this matrix is positive definite in 
the stable region < 2, so an implicit function argument shows that we can define k as a function of ui. 
The second equation is equivalent to the condition n_ (D) — 0, while the third condition is equivalent to to 
the condition that n_ (D) = n_ (— J). Finally the last equation restricts to the subspace (lu, 1) = 0. ■ 
Given this algebraic representation of the stable set it is relatively straightforward to check that the 
region is convex. 

Theorem 2.5 The set is convex. 

Proof. We will first show that the Hessian of m (with respect to u>l) is negative definite for each i. This 
implies that — is a convex function. This implies that t, being a sum of convex functions, is convex. Finally, 
the region = {us \ t(u>) < 2}, is bounded by a level set of a convex function and is therefore a convex 
set. 

Differentiating \2.21\ twice gives 



d m , dm dm 

"5 a !~ 2 Q Q !~ 2 °k,iOk,l 

OLOkduji duj k duji 



N 

E 

3=1 



d 2 *3 
du k duji 



Recalling the definition of M in [ 2.22) , 1 2.23) can be written as a matrix equation in the form 

M(K)h M 

where we have the unknown TV-vector 

hki '■ 



dm dm 

2t> — q 1- 25k,iOkj, 

dbj k dui 



(2.23) 



(2.24) 



d 2 m 



dujkduii 



N 



To clear up a potential confusion: here we think of ( 2.23) as an equation where we fix k, I which determines 
a vector where we vary the functions m in the entries of the vector. 
If M(k) is invertible, then 

h kl = M- 1 ( 2^ ^ + 2S k ,A,) . (2.25) 
\ duj k dui ' J 

We have already shown that M is invertible in the proof of Lemma 2.4 above, and now we compute the 
inverse exactly. Write M = D + 1 (g> 1*, where D is the diagonal matrix whose ith entry is —2m- Using the 
Aronszajn-Krein formula [], or checking directly, we have that 



M = D 



(D _1 1) <g) (D" 1 !.)* 
l + (l,D-il) 



(2.26) 
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In this case, we have 

1 - 1 

((D- 1 l)®(D- 1 l)% = 5r;ri 1 + (1,D- 1 1) = 1 + ^— , (2.27) 



i=l 



so that 

(M-% = -^ 



2«i l + Etia^^/t/ 
In particular, every element of M _1 is negative for k e <S W by < |2.19) . Then we can write 



This is a negative definite matrix: it is given by linear combination of N + 1 matrices with negative co- 
efficients, and each of these matrices are manifestly positive definite — the first is diagonal with positive 
entries, and the rest are of the form (Vk, ) <g> (Vkj). 

Finally, we note that if «j has a negative definite Hessian, then n^ 1 has a positive definite Hessian. To 
see this, write Vi = \j m, and we have 

Since «j > 0, this implies that the Hessian of 2/j is positive definite. Since t = ^ k^ 1 is convex the stable set 
Sg, which is bounded by a level set of r, is convex. ■ 

2.3 Lower Bounds on the stable set 

We next characterize a subregion of Sg for which we are guaranteed stability. The first region takes the form 
of a curvilinear polytope, and is the image of a cube under the map f . We will show that the vertices of this 
polytope actually lie on the boundary of Sg as well. This region is more difficult to deal with analytically, 
so we introduce a second (flat) polytope which is determined by the convex hull of the vertices of the 
curvilinear polytope. This polytope is the Voronoi cell of the root lattice A n and as such its properties are 
well-studied []. This allows us to get very good estimates geometric quantities such as the volume, the 
probability that a random vector lies in the stable synchronous region, etc. 

Definition 2.5 We define the stable cube C to be 

C := {6 \ e [0,71-/2]^} . (2.30) 

and the set V to be the vertex set of the the cube with two vertices removed: 

V := {0|Vi, 9t E {0, ^/2}}/{(0, 0, 0, . . . ,0), (|, |, |, . . . , |)}. (2.31) 

Remark 2.3 One can consider vertices of the n-cube as representing the partition of the n oscillators into two sets. In 
this interpretation V represents the partition of the n oscillators into two non-empty sets. We will see below that the set 
V represents points in T N that map to points on the boundary ofS^. The points (0, 0, 0, . . . , 0) and (§)§)§,..-,§) 
map to the interior ofS^ (in fact to the origin - the most stable configuration). 

Lemma 2.6 The Jacobian is positive semi-definite on the cube C. Moreover, 

dim(kcr(J)) = I*' 6 e ^ C J V > (2.32) 
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Proof. If 8i e (0, it/2), then 6 t - 9j e (-tt/2, tt/2) for all i, j and thus cos(6» l - 0j) > for all i, j. Then 
the Jacobian matrix takes the form of a weighted graph Laplacian on the complete graph with positive 
weights. (Said another way, J has zero row sums and positive off-diagonal entries.) The dimension of the 
kernel of a graph Laplacian is equal to the number of connected components of the graph [] . Since none of 
the weights vanish, the dimension of the kernel is equal to 1. 

When is on the boundary of C, some of the weights may vanish, if the corresponding angles differ by 
exactly tt/2. At the level of the graph, this amounts to decomposing the vertex set of the graph into two sets 
(one corresponding to 0i — and the other to (9; = Z) and breaking all connections between the two sets. 
For the complete graph on N nodes, this will disconnect the graph only if all points belong to one of these 
sets and neither one is empty. These configurations correspond exactly to points in V. ■ 

We have defined a cube in configuration space where the stationary solutions are at least marginally 
stable, and the vertices of this region actually lie on the stability boundary. Next we characterize the shape 
of the corresponding region in frequency space. 

Proposition 2.7 The image ofV under the map f gives 2 N — 2 distinct frequency vectors. These vectors form (up 
to scaling) the vertices ofV(A n ), the Voronoi cell of the root lattice An- Equivalently V(An) is the projection of the 
cube [— onto the (N — l)-dimensional plane normal to 1 in R N . 

Proof. It is easy to check directly that the image of V under the map f consists of vectors of the following 
form: if the vertex has i angles of and j angles of f then the image is a permutation of the vector 

oj = i,i .., i , -j, ~j, —j, -j ) 

j times i times 

where i + j = N. It is a reasonably well-known fact (see, for instance, Conway and Sloane |5| Chapter 
21. 3B) that the Voronoi cell around the origin of the root lattice An has vertices given by vectors of the form 

V l 7V'iV'iV'"''iV' JV" AT' JV"""' N' 

S V ' v ' 

j times i times 

and is the projection of the unit cube \—\, \] N onto the (N — l)-plane of mean zero vectors. 

■ 

Remark 2.4 The root lattice An and the associated polytope arise in a surprising number of areas of mathematics 
including Lie algebras and root systems, Coxeter groups, coding theory, etc. The appearance here is perhaps not 
so surprising since the symmetry group of the Kuramoto problem, Sn X S2 (corresponding to permutation of the 
oscillators and 9 n- —6), is the same as the symmetry group of the A N lattice. 

Corollary 2.8 cl(5 w ) contains V(A n ) and the image of the vertex set f (V) lies on the boundary ofS^. 

Proof. This is clear from the proposition above: a convex polytope is equal to the convex hull of its vertices, 
and the convexity of implies that the convex hull of a collection of boundary points is contained in the 
closure. ■ 
It is worth noting that there is a relatively easily computable set of points lying on the boundary on the 
stability region. While we have not had occasion to use this fact in the present paper this fact may prove 
useful for improving the current bounds, and is presented here. 

Proposition 2.9 The stable region contains the (suitably scaled) dual polytope to the Voronoi cell V(An). The 
vertices of this dual polytope are given by the vectors 

±6^,(1,0, 0,...,-l) 

and all permutations, where lon is the number 

u N := (V32 + (iV-l)2 + 3(7V - 1)) \JlQ + (N — l)y/32 + (N — l) 2 — (N — l) 2 
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Proof. The configurations which correspond to these frequencies are those with a group of N — 2 angles 
at the origin and two others offset symmetrically from these: 

-x, i = 0, 
0, ie(2,N- IS. 
x, i = N 

The corresponding frequency vector is 

uo t = {{N - 1) sin(z) + sin(2ac), 0, 0, . . . , 0, -(N - 1) sin(a;) - sin(2a;)). 
If one maximizes the length of this frequency vector over all x then one finds the expression above. 



3 Upper and Lower Bounds on 7 

We first consider a formal argument as to the sizes of the vertex set in frequency space, i.e. the set f (V). 

Definition 3.1 We define the frequency vectors u> mm , w max as follows: 

« m i„:=±(l,l,...,l,-(JV-l)) f (3.1) 

for all N. We define w max differently if N is even or odd. If N is even, we define u> max as that vector whose first N/2 
components are N/2, and whose last N/2 components are —N/2. If N is odd, we define u; max to be the vector with 
the first (N — l)/2 entries are (N + l)/2, and whose remaining entries are—(N— l)/2. In summary, 



Lemma 3.1 We have 



' (N/2, N/2, N/2, -N/2, -N/2, -N/2), N even, 

v ' " v ' 

N/2 N/2 

(N + l)/2, ...,(N+ l)/2,—(N - l)/2, ...,-(N- l)/2, N odd. 

v ' s v ' 

(JV-l)/2 (N+l)/2 



,- ., , . m i' j /'i i ii— f^ 3 /4, N even, 



(3.2) 



min||f(0)|r=JV(iV-l), max||f(6>)H = { I w ' (3.3) 



eev eev 



lA^(A^ 2 -l)/4, N odd. 



Moreover, the u> which minimize (resp. maximize) in (3.3 1 are permutations ofu) mm (resp. o; max/ ). 

Proof. First note that |V| =2^ — 2 since we choose all vertices of the cube except two. Given 6 £ V, let 
i G 1, . . . , N — 1 be the number of entries of 9 which are equal to and j = N — i be the number which are 
equal to n/2. Everything is the same up to permutation, so replace with the vector 

6* = (0,0,...,0,7r/2,7r/2,...,7r/2) t . 

s * ' y v ' 

i j 

Then 

u>* = f(0*) - -i) * (3.4) 

i j 

We have 

\\u*f =i 2 j + ij 2 =ijN. 
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Clearly, to minimize this, we choose i = N — 1 and j = 1 (or vice-versa), and this gives u>* = u} mm , and 

H^min|| 2 =iV(iV -1). 

To maximize this, it depends on the parity of N; if N is even, we choose i = j = N/2, and if N is odd, 
we choose i = [N/2\ = (N — l)/2. In either case we obtain u max as defined above. We calculate that for TV 
even, ||u; max || 2 = iV 3 /4, and for N odd, ||uw|| 2 = N(N 2 - l)/4. 

■ 

Remark 3.1 Notice that, for all N, ||u; max || 2 = N 3 /4+o(N 3 ). The 0's corresponding to Lj mm ,uj m3X have a relatively 
simple characterization. The minimum vertex can be described as "a herd of sheep and a lone wolf": (N—l) oscillators 
at 9i = and 1 oscillator at Bi = |. The maximum vertex can be described as the state of "two competing cliques" . 
Finally, note that in any case, the vertices are really only distinguished by the size of the partitions i and j, and thus 
there are L^2~^J different types of vertex. 

Lemma 3.2 We have 

1 1 

7min(A0 = mOU\ ' 7max(A0 = T~ f \tf7&T\ > ( 3 - 5 ) 

wpeese ll f (0)ll mfee5 e p(0)|| 

Proof. Recall the definition of r y m \ n (N) , j m3X (N): 



7min (iV) := inf 7*(w), 7max (A0 := sup 7 *(w) g7| 
M=i |w|=i 

Moreover, notice that multiplying the right-hand side of |1.2| by a scalar does not change anything about the 
existence or stability of a fixed point. Choose ||w| = 1, and if we try to solve 

7 f(0) = -u, 

we have to have 7 || f (0)|| = 1, or 

1 

7 " su Pfle58 ||f(e)ir 

Since Sg is open, this bound is saturated and gives us 7m in(^V)- The opposite argument works for 7max (-/V). 
■ 

Theorem 3.3 Recall the definition of 7min , 7max in \1.7) above. Then 

(2N- 3 / 2 , N even, 

7min( >~\2N-V 2 + 0(N- 5 / 2 ), Nodd, { ' 

and 

1 /2 
; < 7max (iV) < , (3.7) 

^N(N-1) " 7maxl y/N(N-l) 

Proof. We start with j m \ n (N). Recall the identity ( |2.21^ ; summing both sides of this equation over i 
gives us 

N N N N 

i—l 2—1 i— 1 i — 1 

The right-hand side of \3.8\ is maximized when we choose Kj = N/2 and therefore we have the global 



bound ||w| < 7V 3 /4. From Lemma 3.2 we can deduce that ~/ m \ n (N) > 2/N 



-3/2 



If N is even, then this bound is obtained at o; max , and therefore 



9 iV 3 

sup IMI 2 = — , (3.9) 

u6S„ 4 
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and thus 7 min (iV) = 27V- 3 / 2 . 

If N is odd, then uj max does not attain the N 3 /A bound. However, because of o; max , we know that 

9 N 3 - N 
sup ||f(0)|| 2 > — — , 
ees e 4 

and therefore 

7min (iV) < = ^^j^ = ^(1 + O(N-i)). 

In a similar vein, we have that 

inf ||f(0)|| < min||f(0)|| = y/N(N-l), 

"Gog fGV 

and so 7 max (N) > l/(y/N(N- 1)). The upper bound comes from the fact that the stability region S^, 
contains the Voronoi polytope. It is straightforward to calculate the radius of the sphere inscribed in the 
Voronoi polytope: it is ^ independent of dimension n. This gives a lower bound for infg f (w), and thus an 
upper bound on -f max of 

7m ax(A0 < ^^TJ 



Remark 3.2 PV/rife £/ze above estimates give the correct order of magnitude, and in certain cases the exact value, of 
7 m j n and 7 max we believe that the following are the exact values. 



7max(A0 = 



1 



V N ( N - 1) 

2 

AT 2 



7min(#) = <| 8^/2(^-1) 

(V8iV 2 - 16JV + 9 + 3) v^iV 2 - 8iV + 3 + vW 2 "- 16iV + 9' 



iV even, 
N odd. 



We conjecture that the configuration with N — 1 oscillators having angle 6i — and one having angle 9i = | 
is a global minimizer of the frequency w over the marginally stable set. It is easy to check that it is a local minimum, 
and numerical results indicate that (at least for small N) it is a global minimum. For j m \ n the even case is tight, 
from the calculation above. The conjectured value for N odd requires some comment. For the case N odd there is a 
configuration which is not a vertex which has a larger value ofuj than any vertex. The configuration which gives this 
is as follows: a single oscillator at 6 = 0, a two groups of ^=1 placed symmetrically on either side at angle x. This 
gives the magnitude of the frequency as the solution to the following maximization problem 

N - 1 



lo* = sup y/N - l(sin(x) H — sin(2ar)) 

X ^ 

whose solution is 



(V8N 2 - 16N + 9 + 3) V^N 2 - SN + 3 + VM 2 - 16N + 9 
8y/2(N-l) ' 

3 

It is straightforward to check that for large N this is asymptotically 4^- + 0{N^), and would give j m m(N) 



2N 3 / 2 + 0(N r °/ 2 ), consistent with Theorem 3.3 
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4 Large N limit 

In many of the problems of interest one is interested in the case where the number of oscillators is large, and 
the frequencies are chosen from a specific probability distribution. In this section we establish rigorously 
that the interesting scaling for the Kuramoto problem when the frequencies are chosen independently is 

not the classical scaling TV -1 , but actually the scaling <p(N) := ^/V+i^ . More specifically, we prove the 
following: 

Theorem 4.1 Suppose that the frequencies are independent identically distributed (i.i.d.) Gaussian random vari- 
ables with unit variance, i.e. we assume that the LOi are chosen independently, and that 

P(wi G (0,6)) = — L / e- x2 ' 2 dx. (4.1) 

V2-7T J a 

Let Vsyndl, N) denote the probability that the system \\.2) has a stable state with all oscillators locked. Then we 
have the following dichotomy: 

• If 5 < 1 then lim^oo V sync {5(p(N),N) = 0. 

• If 8 > 2 then limjv-,00 V sync (6ip(N) , N) = 1. 

The basic strategy is straightforward: given our previous results on the shape of the stability domain 
we establish upper and lower bounds on the probability that a Gaussian random frequency will lie in the 
stable region. We prove each of these statements separately in two lemmas. 



Lemma 4.2 Suppose that the components of the frequency are i.i.d Gaussian as in | |4.1| . The probability that the 
system will exhibit stable synchronization satisfies the upper bound 

Psynd^N) <^ (erf (7JV/V2))* \ (4.2) 
Proof. The vector u> is distributed according to the multivariate Gaussian measure 

P(w G A) = (2tt)-^ / exp(- |x| 2 /2) dx. (4.3) 

J A 

Since we have moved to the co-rotating frame and we have assumed that J^^i = 0' we ma ke the 
following (non-orthogonal!) change of variables. First choose Xi as m ( |4.3) above, then write 

1 N 

Ui=Xi~v, v: =jj^2xi, i = l,...,N. (4.4) 

Note then that 



N 

i=l 



JV-1 JV-1 

UN = - ^ u>i, or, xn = v - ^ u i- ( 4 - 5 ) 

i=l i=l 

Note that we've written the %'s in terms of uii with i = 1, . . . , N — 1 and v; these will be the new variables. 



As we prove in Lemma 4.3 below, the Jacobian of this change of variables is N. Moreover, the quadratic 
form transforms as 



JV JV-1 



= 1 i=l 

JV-1 JV-1 / N-l 



i=l 

N-l N-l 

J2 <4 + 2« + ( N - ^ y2 + ( v - 12 ^ ) < 4 - 6 ) 

i=l i=l \ i=l 

JV-1 /JV-1 \ 2 

, <=1 
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Therefore, for any set A, if we denote the transformation in 1 4.4 4.5 1 as q, then we have 



-|x| 2 /2 



dx 



l{A) 



Ndu>, 



(4.7) 



where the N comes from the Jacobian of the transformation. To determine the domain of integration A, 
note that if we have a fixed point, 



N = 7 



N 



0=1 



so we take our domain to be 



A := \oji\ < jN for i = 1, .... TV - 1, «eK. 



(4.8) 



(4.9) 



and it is a necessary condition for synchronization that x <= A. Thus we compute an upper bound on the 
probability for synchronization: 



?>sync(7,A0 < (27r)- JV / 2 



N 



[--/N,-yN] N ~ 1 xl 



N Y[ dcJi dv 



i=l 



jV 



< (2n)~ N / 2 N / e-MES'^W] TT du> l dv 

J[-jN,fN] N - 1 xM ., 



N-l / , - 7 iV 

^TT M= / 



e- w '/ 2 dui]x[ 

-jN J \ V^TT J-c 



-Nv 



dv 



N-l 



N Yl erf(77V/v / 2) x 1 = \/Ar (erf^N/V?) 



i=l 



N-l 



Lemma 4.3 If we define the transformation as in {4.4 4.5*, then its Jacobian is 

d(xi, ■ ■ -,Xn) 



Proof. We compute 



dXr 



diOi dv 



5(wi, . . .,u N -i,v) 



= N. 



1, t,j = l,...,JV-l, 



u 3 

9xn 



= 1, ^=1. 



(4.10) 



duii dv 

Thus this Jacobian (which we will call Gjv) is the constant matrix which we write in block-diagonal form: 



Gn — 



In-i 1(jv-i)> 



-1 



lx(JV-l) 



1 



(4.11) 



where Ik is the k x k identity matrix, l aX 6 is the a x b matrix of all ones, etc. By elementary row operations 
(adding the sum of the first N — 1 rows to the last) the Jacobian can be reduced to 

In-i l(JV-l)xl 



Olx(JV-l) 



N 
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This matrix is upper-triangular and we can read off the determinant as N. ■ 
The next proposition gives a lower bound for the probability by identifying a region where a sufficient 
condition for stability holds, over which the Gaussian integral can be evaluated rather explicitly. The proof 
of this is straightforward but requires some facts about polytopes in M. N 



Lemma 4.4 The probability of synchronization satisfies the following lower bound: 



N 



(4.12) 



Proof. The proof here uses the fact that the polytope V(An), which we know to be contained in the 
(closure of) the stable region, is the projection onto the mean zero hyper-plane of the cube in one higher 
dimension. This, together with the orthogonal invariance of the Gaussian, give the result. 

We use the notation of the previous lemma, that x € a frequency vector which is not assumed to 
have zero mean, and u; = x — %^ is the orthogonal projection onto the mean zero subspace. By the 
orthogonal invariance of the Gaussian we have that 



dx 



(27T) 



2 du). 



IV(A N )xR JV(A N ) 

Since the stability region contains the polytope V(A n ), we have the inequality 



Psync(7,A0> / (27T)- 

Jv(A N )xM. 



dx= / (2tt) 

<V(A N ) 



2 dui. 



N 



Next note that since the polytope V(A N ) is the projection of the cube — ^ , ^ onto the (N — 1)- dimen- 
sional plane normal to 1, we necessarily have that the cube is contained in V(A N ) x M, the cylinder in R N 
with cross-section V(Ajy): 

1 V 

2' 2. 



C V(A N ) x 



This in turn shows that 



Vsync(l,N) > 



V(A N )xl 



. . _ jy_ ii jc ii 

(2ir) 2 e 2 dx > 



(2tt) 2 e 2 dx = (erf ( — j=)" . 



Proof of Theorem 4.1 Assume that 7 = 5tp(N), with S < 1. Using Lemma 4.2 we have that 

^sync(7, N) < \fN (erf( 7 iV/\/2)) " ' . 
The standard asymptotic expansion for the error function gives [], for large x: 



14.21 



erf(x) = 1 



;(1 + 0(X- 2 )), 



(4.13) 



so we have 



vi( 7 N/V2) = eri(6VHNj) = 1 - ^^7=^ (1 + ^(^(TV)- 1 )). 



This means that there exist C\ € (0, 1), S < 8 < 1 such that, for N sufficiently large, 



(4.14) 



erf (7AT/V2) = 1 - C 2 exp(-<T ln(N)) < 1 - C 2 iV" 
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A straightforward application of L'Hopital's Rule shows that, if 5 < 1, then 



1 - C 2 N- S ) x exp(-C 3 iV 1 - 5 ), 



and therefore, if 5 < 1, in the limit as N — > oo, the right-hand side of | |4.2) goes to zero as iV — > oo. (The v 
in front and one power of erf in the back will not change anything for large N.) 



Now assume 7 = Sip(N) with S > 2, and the argument is similar. Using Lemma 4.4 we have 



Again using [ 4.13) , we plug 7 = 5tp(N) into a single term and obtain 



f( 7 7V/(272)) = crf^v/L^V)^) = 1 - 6X f ( * lnW) (1 + OiHN)' 1 )) 



or 



vViogfjv) 



This means that there exists C > 0,2 < S < S such that for N sufficiently large, 

erf( 7 7V/(2V2)) > 1 - CN~ p/i . 
Again using l'Hopital's Rule, we have that if 6 > 2, then this goes to 1 as N — s- 00. 

■ 

Finally, we prove that, in the standard scaling, the probability of choosing a fully synchronous solution 
is zero and, moreover, that it goes to zero exponentially fast. Specifically, we prove: 

Proposition 4.5 For all S > 0, there exist C%, C2 eM with 

Vsync(S/N, N) < C ie - C2{5)N . (4.15) 



Proof. Using Lemma 4.2 we have that 

V sync (6/N,N) < ViVerf^/v^ 1 

Choose erf (S/^/2) < a < 1, and write /3 = a/ erf(<5/v / 2). Then a calculus argument shows that if we define 
C*i > (2elog(l//3))- 1 / 2 ,thenCi > v^V^" 1 for all N, or 



erf (6/V2) N < da 1 *- 1 . 
Choose C 2 — — log a and pull a power of a into C\, and we are done. 

■ 

Finally we close this section with some numerical simulations. The figure [T] depicts a Monte-Carlo 
simulation of the Kuramoto problem with N = 1000, N = 2000, N = 10, 000 and N = 25, 000 oscillators. 
The figures were generated as follows: for each realization a direction was generated uniformly on the 
TV — 1-dimensional sphere, and the distance to the boundary of the stability region was computed. In this 
situation the formulation of Mirollo and Strogatz was found to be more computationally efficient, and this 
was solved numerically using a bisection method. For each such direction, given the distance to the stability 
boundary the conditional probability that a Gaussian random vector in that direction would lie within the 
stability region could be calculated, and averaging over all realizations gives a numerical approximation to 
the probability of synchronization. For each of the graphs 5000 realizations were used. 

The numerics agree well with the analytical results. We have shown analytically that in the limit of large 
N the probability of full synchrony is for 7 < 1. The numerics suggest that 

with a critical coupling constant of roughly 7* sw 1. This statement is consistent with, but sharper than, the 
conclusion of Theorem 14. II 
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Figure 1 : A numerical simulation of the probability of full synchrony in the Kuramoto model in the newly 
identified scaling. 

5 Examples 

5.1 Comprehensive example for three oscillators 

For N = 3 there are 2 3 — 2 = 6 vertices. The corresponding frequency vectors are given by 

tfi - (1,1,-2)* 
w 2 = (1,-2,1)* 
c3 3 = (-2,l,l)* 
tf 4 = (-1,-1,2)' 
lu 5 = (-1,2,-1)* 
06 = (2, -1,-1)*. 



The plane orthogonal to the vector (1, 1, 1) is spanned by the vectors 

ei = (l,0,-l)/V2, e 2 = (1, -2, 1)/V6. 

In this basis the frequency vectors have the representations ±\/6e2, \/6(± 2 e * 5 e2 )' which are obviously 
the vertices of a regular hexagon of side length \/6. As shown in the lemma, these represent local minima 
of distance on the surface of marginal stability. 

The phase diagram for three oscillators is summarized in Figure|3] Since the map from the configuration 
space T' i_1 to the frequency space has degree zero it follows that the number of preimages of a given 
point to (the number of configurations with a given frequency) is even, half of which have positive (reduced) 
Jacobian determinant and half of which have negative Jacobian determinant. The latter are, of course, 
always unstable. The former have even index but may not be stable. Outside the shaded region there are 
no fully synchronized solutions. As one crosses the stability boundary a pair of synchronized solutions are 
created: one is stable (index zero) and one is unstable (index one). The guaranteed stable region, the image 
of the cube [0, under the frequency map, is also plotted: it is a dark curve just inside the stable region. 

As shown earlier we can derive expressions for the frequency vectors which are hardest and easiest to 
synchronize. The hardest vectors to synchronize have (n — 1) oscillators traveling together with no phase 
shift and 1 oscillator leading (or trailing) by phase | . These solutions have frequency vectors 

«=(±1,±1,t2) 
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and permutations. These points are marked by the six dots on the boundary. These points lie on a circle of 
radius V&. The solutions which are easiest to synchronize have one oscillator at phase 0, one leading by x 
and one trailing by x with corresponding frequency vector 

u = (±(sin(a;) + sin(2z)), 0, =F(sin(a:) + sin(2a;))) 

plus permutations. The maximum length frequency vector of this form is attained at 

x = arccos(^(\/33 — 1)) 
8 

which gives a frequency vector of length 




It is interesting that, in the case of three oscillators the stability region is very close to a circle and the 
phase transition to the fully-synchronized state is quite sharp: there are no fully synchronized solutions 

for uj > (3 + \/33)^/| (15 + \/33) w 2.49 and all solutions synchronize for co < \f& w 2.45, a relative 
frequency shift of about 1.6%. 

On the interior of the stability region there are two more secondary bifurcation curves, which look like 
a pair of curvilinear triangles rotated through angle ^ relative to one another. As one crosses these curves 
a pair of unstable solutions are created, one of index 1 and one of index 2. Thus the frequency plane can be 
subdivided as follows: 

• Region 1: Two solutions, one of index and one of index 1. 

• Region 2: Four solutions, one of index 0, two of index 1, one of index 2. 

• Region 3: Six solutions, one of index 0, three of index 1, two of index 2. 

Also shown is a related figure in the the configuration space. Here the configuration space is colored 
according to the stability of the given configuration (note that by the 0(1) symmetry the first oscillator can 
be chosen at 6\ = 0). The light-colored central region represents the stable solutions. It is, as was shown, 
convex. The surrounding gray region represents the index 1 solutions, i.e. those with a single unstable 
direction. Both of these regions cover the range of the frequency map. Finally the two darkest regions 
represent the two triangular regions where there exist solutions of index two. 

5.2 Example for four oscillators 

The case of four oscillators is the first with different types of vertices. The vertices are of two types: the first 
are vertices of the form 

oJ=(±1,±1,±1,t3) 

and permutations, which have length ^12. These are local minima of the length of the frequency vector 
constrained to the surface of marginal stability, and represent the hardest frequencies to stabilize. There are 
eight such vertices. The second type of vertex has frequency vectors of the form 

lu = (±2,±2,T2,T2) 

and permutations, which have length 4. There are six such vertices and these represent the easiest vectors 
to stabilize. These points together form the vertices of a polytope known as the rhombic dodecahedron. It 
is the Voronoi cell for the ^ lattice - the face centered cubic lattice. 

The figure depicts the region of guaranteed stability - the image of the cube [0, f] under the frequency 
projection map. Since this is a symmetric nonlinear projection of the cube in K 4 into R 3 it is perhaps not 
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Figure 2: A plot of the stability regions for the three particle system. In the light beige central region the 
Jacobian is of index 0. In the surrounding pale blue region the Jacobian is of index 1. In the two dark 
blue islands the Jacobian is of index 2. The marked points depict special points. The six central magenta 
ones depict the last frequencies to be stabilized, while the six green ones represent points on the boundary 
between the regions of index 2 and index three. 

surprising that it takes the form of a curvilinear rhombic dodecahedron, since the corresponding linear pro- 
jection gives the (flat) rhombic dodecahedron. The vertices typified by the frequency oj = (1,1,1,-3) are 
those where three of the curvilinear rhombic faces come together, while those vertices typified by frequen- 
cies like (2, 2, —2, —2) are those where four rhombic faces meet. The actual stability region (not depicted) is 
somewhat larger, and resembles an octahedron which is tangent to the region of guaranteed stability at the 
fourteen vertices listed above. 



In the preceeding section we showed that the critical scaling for full synchronization differs from the clas- 
sical scaling by a logarithmic factor: that one should consider 



the correct scaling, which is connected with the extreme value statistics of the frequencies |4}|23). 

Let us recall the basics of extreme value statistics. Give a collection of idenpendent and identically 
distributed random variables u)i the extreme values statistics concerns the distribution of the quantity Mat = 
max(cji, u>2, W3, . . .u>n). The Fisher-Tippett-Gnedenko theorem J23) characterizes the rescaled distribution 
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Figure 3: The frequency map. The magenta points correspond to the analogous points in the previous plot: 
the boundary of the region between the region of index zero and the region of index one corresponds to 
the boundary of the range. These points denote the last frequency vectors to exhibit synchrony. The green 
points denote similarly distinguished points on the boundary between the regions of index 1 and index 2. 
Also shown is the guaranteed stable region, which is tangent to the boundary at the magenta. 

of such a quantity. It says that if the distribution of a suitably rescaled Mjv converges to a non-degenerate 
distribution G(z): 

P(Mjv - b N )/a N < z)^G(z) as TV — > oo 
then G{z) is one of the following three distributions, the Gumbell, Frechet and Weibull distributions. 

Gi{z) = exp(exp(-z)), (6.1) 
nM fo, z < 0, 

G2(z) = z~>o, (6 - 2) 

Gs(z) = h P( -^' Z< ° n (6.3) 
II, z > 0. 

We claim that the preceeding calculation shows that the probability of full synchrony for the Kuramoto 
models is bounded above and below by extreme value statistics for the Gaussian, and that the slightly 
anomolous scaling seen in this problem is a reflection of the extreme value statistics. 

It is easy to see that for the Kuramoto problem one has the following obvious estimate for synchronous 
solutions. From the formula for the frequencies in the classical scaling one has 

uii - uij = — ^2 sm(9 k —0i)- sin(6> fc - 9j) 

k 
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Figure 4: A region of guaranteed stability for the Kuramoto problem with four oscillators. The region is a 
rhombic dodecahedron with curvilinear faces. The vertices are all points of marginal stability, where there 
is a zero eigenvalue of multiplicity two or more. 

leading to the easy estimate 

max(wi) — min(wi) < 2 

i i 

that must hold in order for there to be a fully synchronous solution. This gives an upper bound for the 
probability in terms of the distribution of 

max(u)J — min(a;i), 

i i 

a kind of mean adjusted extreme value statistic. This holds for any choice of distribution of the frequencies. 
For the Gaussian case it is clear that the correct scaling is for this statistic is 

max(wj) — min(o;j) oc ylog N. (6.4) 

i i 

This gives a motivation for the modified scaling necessary for full synchrony. For Gaussian distibutions of 
frequencies the lower bound is also of the same form. The lower bound we proved in the preceeding section 
estimates the probability of synchrony in terms of the probability that the frequency lies in a Voronoi cell of 
the Am lattice. 

So then might pose the question of where the standard scaling might be applicable? We note that the 
anomalous scaling occurs because of the fact that while the typical sample of a Gaussian is O(l), the max- 
imum of a sample of size N (the first order statistic) is significantly larger than this, being 0(^/log N). So 
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one could modify the original question: instead of requiring full synchrony for all oscillators, which leads 
to a condition like l |6.4| , we could ask the question of whether or not we have "all but one" oscillator syn- 
chronize, or even "all but k" oscillators synchronize. This would then be governed by the typical size of 
the second order statistic or the k + 1st order statistic. However, for fixed k and TV — > oo, these all have the 
same scaling, so we conjecture that for these problems we would also obtain the scaling seen here. 

However, if we modified even further and asked for solutions such that some fixed fraction of oscillators 
synchronized (e.g. for large N, we require that 0.97V oscillators synchronized), then this would come from 
the typical size of the 90th percentile of the sample, and this is 0(1). Thus we further conjecture that this 
type of question would lead to the classical scaling. We will consider these questions in future work. 

Acknowledgments 

The authors would like to thank Yulij Baryshnikov for useful discussions. RELD was partially supported 
by NSF grant CMG-0934491. JCB and MJP were partially supported by NSF grant DMS-0807584. 

References 

[1] J. A. Acebron, L.L. Bonilla, C.J.P. Vicente, F. Ritort, and R. Spigler. The Kuramoto model: A simple 
paradigm for synchronization phenomena. Reviews of modern physics, 77(1):137, 2005. 

[2] Neil J. Balmforth and Roberto Sassi. A shocking display of synchrony Phys. D, 143(l-4):21-55, 2000. 
Bifurcations, patterns and symmetry. 

[3] Matthew Bennett, Michael F. Schatz, Heidi Rockwood, and Kurt Wiesenfeld. Huygens's clocks. R. Soc. 
Lond. Proc. Ser. A Math. Phys. Eng. Sci., 458(2019):563-579, 2002. 

[4] Eric Bertin and Maxime Clusel. Generalized extreme value statistics and sum of correlated variables. 
/. Phys. A, 39(24):7607-7619, 2006. 

[5] J. H. Conway and N. J. A. Sloane. Sphere packings, lattices and groups, volume 290 of Grundlehren der 
Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences]. Springer- Verlag, New 
York, third edition, 1999. With additional contributions by E. Bannai, R. E. Borcherds, J. Leech, S. P. 
Norton, A. M. Odlyzko, R. A. Parker, L. Queen and B. B. Venkov. 

[6] G. Bard Ermentrout. Synchronization in a pool of mutually coupled oscillators with random frequen- 
cies. /. Math. Biol, 22(l):l-9, 1985. 

[7] E. Fermi, J. Pasta, and S. Ulam. Studies of nonlinear problems. Los Alamos document LA 1940, 1955. 

[8] J. Ford. The Fermi-Pasta-Ulam problem: Paradox turns discovery. Physics Reports, 213(5):271-310, 
1992. 

[9] Seung-Yeal Ha, Eunhee Jeong, and Moon-Jin Kang. Emergent behaviour of a generalized Viscek-type 
flocking model. Nonlinearity, 23(12):3139-3156, 2010. 

[10] Seung-Yeal Ha, Corrado Lattanzio, Bruno Rubino, and Marshall Slemrod. Flocking and synchroniza- 
tion of particle models. Quart. Appl. Math., 69(1):91-103, 2011. 

[11] D. Hansel and H. Sompolinsky. Synchronization and computation in a chaotic neural network. Phys. 
Rev. Lett., 68(5):718-721, Feb 1992. 

[12] C. Huygens. Horoloquium Oscilatorium. Parisiis, Paris, 1673. 

[13] B. W. Knight. Dynamics of encoding in a population of neurons. Journal of General Physiology, 59(6):734- 
766, 1972. 



24 



[14] N. Kopell and G. B. Ermentrout. Symmetry and phaselocking in chains of weakly coupled oscillators. 

Comm. Pure Appl. Math., 39(5):623-660, 1986. 

[15] Y. Kuramoto. Chemical oscillations, waves, and turbulence, volume 19 of Springer Series in Synergetics. 
Springer- Ver lag, Berlin, 1984. 

[16] Y. Kuramoto. Collective synchronization of pulse-coupled oscillators and excitable units. Physica D, 
50(l):15-30, May 1991. 

[17] Georgi S. Medvedev and Nancy Kopell. Synchronization and transient dynamics in the chains of 
electrically coupled FitzHugh-Nagumo oscillators. SIAM J. Appl. Math., 61(5):1762-1801 (electronic), 
2001. 

[18] Renato E. Mirollo and Steven H. Strogatz. Synchronization of pulse-coupled biological oscillators. 
SIAM }. Appl. Math., 50(6):1645-1662, 1990. 

[19] Renato E. Mirollo and Steven H. Strogatz. The spectrum of the locked state for the Kuramoto model 
of coupled oscillators. Phys. D, 205(l-4):249-266, 2005. 

[20] C. S. Peskin. Mathematical aspects of heart physiology. Courant Institute of Mathematical Sciences New 
York University, New York, 1975. Notes based on a course given at New York University during the 
year 1973/74, see http : / /math . nyu .edu/faculty/pe skin/ heart notes/ index . html. 

[21] A. Pikovsky M. Rosenblum, and J. Kurths. Synchronization: A Universal Concept in Nonlinear Sciences. 
Cambridge University Press, 2003. 

[22] Michael Reed and Barry Simon. Methods of modern mathematical physics. IV. Analysis of operators. Aca- 
demic Press [Harcourt Brace Jovanovich Publishers], New York, 1978. 

[23] Sidney I. Resnick. Extreme values, regular variation and point processes. Springer Series in Operations 
Research and Financial Engineering. Springer, New York, 2008. Reprint of the 1987 original. 

[24] Barry Simon. Spectral analysis of rank one perturbations and applications, 1993. Lectures at the 
Vancouver Summer School in Mathematical Physics. 

[25] Steven Strogatz. Sync: The Emerging Science of Spontaneous Order. Hyperion, 2003. 

[26] Steven H. Strogatz. From Kuramoto to Crawford: exploring the onset of synchronization in popula- 
tions of coupled oscillators. Phys. D, 143(l-4):l-20, 2000. Bifurcations, patterns and symmetry. 

[27] Dane Taylor, Edward Ott, and Juan G. Restrepo. Spontaneous synchronization of coupled oscillator 
systems with frequency adaptation. Phys. Rev. E (3), 81(4):046214, 8, 2010. 

[28] Mark Verwoerd and Oliver Mason. Global phase-locking in finite populations of phase-coupled oscil- 
lators. SIAM }. Appl. Dyn. Syst., 7(1):134-160, 2008. 

[29] Mark Verwoerd and Oliver Mason. On computing the critical coupling coefficient for the Kuramoto 
model on a complete bipartite graph. SIAM }. Appl. Dyn. Syst., 8(l):417-453, 2009. 

[30] Arthur T. Winfree. The geometry of biological time, volume 12 of Interdisciplinary Applied Mathematics. 
Springer- Ver lag, New York, second edition, 2001. 



25 



